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We show that the acceptance probability for swaps in the parallel tempering Monte Carlo method 
for classical canonical systems is given by a universal function that depends on the average statistical 
fluctuations of the potential and on the ratio of the temperatures. The law, called the incomplete 
beta function law, is valid in the limit that the two temperatures involved in swaps are close to one 
another. An empirical version of the law, which involves the heat capacity of the system, is developed 
and tested on a Lennard- Jones cluster. We argue that the best initial guess for the distribution of 
intermediate temperatures for parallel tempering is a geometric progression and we also propose a 
technique for the computation of optimal temperature schedules. Finally, we demonstrate that the 
swap efficiency of the parallel tempering method for condensed-phase systems decreases naturally 
to zero at least as fast as the inverse square root of the dimensionality of the physical system. 
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Keywords: Monte Carlo, parallel tempering, replica exchange, acceptance probability 



I. INTRODUCTION 

Ergodic Monte Carlo simulations of large dimen- 
sional systems having complicated topologies with 
many disconnected local minima are difficult com- 
putational tasks, though indispensable for many 
physical applications, 1 . 2 . 3 . 4 . 5 . 6 . 7 . 8 ! 9 : 10 . 11 ! 12 ! 13 ! 14 : 15 ! 16 ! 17 . 18 
Among the various methods dealing with such problems, 



a derealization parameter, as in the q-jumping Monte 
Carlo method^ 7 , or suitable modifications of the poten- 
tial, as in the Hamiltonian replica exchange methodiiSii 8 " 
In parallel tempering, swaps involving two tempera- 
tures (3 and (3 1 are attempted from time to time in a 
cyclic or random fashion and accepted with the condi- 
tional probability 



the parallel tempering methoc 



one of the most 



successful, especially given the simplicity of the idea and 
the ease of implementation. For sure, the idea of cou- 
pling two independent Markov chains characterized by 
different parameters in order to ensure transfer of infor- 
mation from one to the other has a long history. In phys- 
ical sciences, coupling strategies have been employed for 
the development of specialized sampling techniques such 
as replica Monte Carlo sampling of spin glasses, 1 jump- 
walkingjSi* 2 ^ and simulated tempering)2i24 to give a few 
examples. How this coupling must be performed in the 
simplest, most general, and most efficient way is, how- 
ever, a quite difficult problem. 

The parallel tempering method, as we utilize it in this 
article, addresses the question of coupling independent 
Monte Carlo chains that sample classical Boltzmann dis- 
tributions for different temperatures and which are usu- 
ally generated by the Metropolis et al algorithmiS 5 ^ The 
method has been formalized seemingly independently by 
Geyer and Thompson 1 ^ as well as by Hukushima and 
Nemoto^S Of course, it is not necessary that the distribu- 
tions involved in swaps differ through their temperature. 
For instance, the controlling parameter may be the chem- 
ical potential, as in the hyperparallel tempering method, 6 



lin |l, ( 



,(/3'-/3)[V(x')-y(x)] 



(1) 



where V(x) is the potential of the physical system un- 
der consideration. This acceptance rule ensures that the 
detailed balance condition 2 ^ is satisfied and that the equi- 
librium distributions are the Boltzmann distributions. 
As Eq. Q suggests, the efficiency of the temperature 
swaps depends on the difference between the inverse tem- 
peratures /3 and f3' . In order to maintain high acceptance 
ratios, only swaps between neighboring temperatures in 
a given schedule j3 min — fli < /3 2 < • • • < (3 N — P max are 
attempted. An optimal schedule of temperatures has the 
property that the acceptance ratios between neighboring 
temperatures are equal to some predetermined value p, 
value that is usually greater or equal to 0.5. The deter- 
mination of the optimal schedule is complicated by the 
fact that the distributions of the coordinates x and x' 
are also temperature dependent (they are, of course, the 
Boltzmann distributions at the temperatures (3 and /?', 
respectively) . 

In this article, we attempt to answer the following im- 
portant question: What are the main properties of the 
physical system that control the acceptance ratio in the 
limit that the difference (3' — (3 is small? The answer 



2 



to this question allows for a better understanding of the 
applicability as well as the limitations of the parallel tem- 
pering method. In addition, it allows for the development 
of optimal temperature schedules in a way that seems 
more direct and easier to implement than other adaptive 
strategies^ 

In Section II, we demonstrate in a rigorous mathemat- 
ical fashion that the acceptance probabilities for parallel 
tempering swaps are controlled, within an 0{\(3' — /3| 3 ) 
error, by the ratio of the two temperatures involved in 
swaps and by the average potential fluctuations of the 
system, at the inverse temperature (3 — {[3 + f3')/2. The 
acceptance probabilities are well approximated by the so- 
called incomplete beta function law, which has the addi- 
tional property that it is exact for harmonic oscillators. 
Under the assumption that the relation between the av- 
erage fluctuations and the average square fluctuations of 
the potential is roughly the one for harmonic oscillators, 
we develop an empirical version of the incomplete beta 
function law, version that connects the acceptance proba- 
bilities for parallel tempering swaps with the temperature 
ratios and the heat capacity of the system. 

In Section III.B, we show how the incomplete beta 
function laws can be employed for the determination of 
optimal temperature schedules. We also explain the em- 
pirical observation that a geometric progression is the 
best schedule for systems and ranges of temperatures 
for which the heat capacity is almost constant^ In Sec- 
tion III.C, we demonstrate rigorously that the efficiency 
of the parallel tempering method for harmonic oscilla- 
tors decreases naturally to zero at least as fast as the 
inverse square root of the dimensionality of the physical 
system. We then argue that the loss in efficiency is even 
greater for condensed-phase systems (both solids and liq- 
uids). This result seems to be in contradiction with the 
findings of Koike^ who suggested that such a curse of 
dimension does not appear for parallel tempering. How- 
ever, the result is in agreement with the explanation of 
Fukunishi, Watanabe, and Takada^ It is for this reason 
that we insist that our findings be proven in a mathe- 
matically rigorous way. 

The rigorous form of the incomplete beta function law 
involves the average potential fluctuations at certain tem- 
peratures. An evaluation of this property by Monte Carlo 
simulations would require the computation of a double 
integral over the configuration space. For this reason, 
current Monte Carlo codes would have to be modified 
extensively in order to take advantage of the incomplete 
beta function law for the design of optimal temperature 
schedules. To circumvent this undesirable situation, we 
propose an empirical version of the incomplete beta func- 
tion law, version that is derived under the assumption 
that the relation between the average fluctuations and 
the average square fluctuations of the potential is roughly 
the one for harmonic oscillators. In Section IV, we illus- 
trate the good applicability of the empirical law by per- 
forming a Monte Carlo simulation for a cluster made up 
of 13 atoms of neon that interact through Lennard- Jones 



potentials. 



II. THE INCOMPLETE BETA FUNCTION LAW 

Consider a d-dimensional physical system described by 
the potential V(x), which is assumed to be bounded from 
below. To simplify notation, we may also assume that 
the global minimum of the potential is 0, perhaps after 
addition of a constant. Clearly, the addition of a constant 
does not change the acceptance probabilities for swaps. 

In the parallel tempering algorithm, swaps involving 
two temperatures (3 and (3' occur with the conditional 
probability 



The joint probability distribution density of the points x 
and x' in an equilibrated system is given by the formula 



1 



Q(J3)QW) 



e -0V(x) e -0'V(x.') 



where 



QW) 



is the configuration integral. It follows that the accep- 
tance probability Ac(f3, (3') for swaps between neighbor- 
ing temperatures is given by the average 



1 



x min < 1, e 



d*j dx'e-^ x V' 3 ' y ( x '> 

|l je (/3'-/3)[V(x')-V(x)]| _ ( 2 ) 



By construction, Ac{(3, (3') is symmetrical under ex- 
change of variables. Without loss of generality, we may 
assume j3' > j3. Then, 



min{l,e (/3 '- /3)[y(x ' ) - l/(x)1 } = e (P'-P)™™{0,V(*')-V(x)} 

= e V^[V{*')-V{K)] e - { -S^l\V{*')-V(K)\^ 

Replacing Eq. J3J in Eq. (J2J, we obtain 



(3) 



Ac((3,(3') = 



1 



Q(P)Q(J3') 

x exp 



dx / dx'e-^ x )e-^ x ') 
R-l 



R+l 



0\V(x')-V(x)\ 



(4) 



In Eq. @, the variables 7? and (3 are defined by the equa- 
tions 



~f3 — 2 — : 



(5) 



respectively. Due to the nature of the results we prove, 
it is more convenient to express the various formulas in 
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terms of the new variables R and (3. Because j3' > (3, we 
need only consider the case R > 1. 

As announced in the introduction, we are interested in 
establishing asymptotic laws in the limit that R — > 1 for 
which the error is of order O (\/3' — /3| 3 ) or, alternatively, 
O (\R — 1| 3 ) . At this point, let us see that the acceptance 
probability given by Eq. Q) is alternatively given by the 
formula 



Ac{p,p') = 



exp 



R- 1 
R+ 1 



B\u'-u\ 







where, in general, (f(U,U'))g denotes the statistical av- 
erage 



1 



CO rOO 



e-^ u e- 0U 'il{U)fl(U')f(U, U')dUdU'. 

(7) 

In Eq. J7J), fl(U) denotes the density of states. 

In these conditions, the following proposition holds 
true. 



Proposition 1 We have 
R-l 



where 



1 



R 



jM (p) + 0(\R~ 1| 3 ) , (8) 



M0)=~H\U' -U\)- p . 



(9) 



Proof. A Taylor expansion of the exponential function 
to the third order and the identity (U 1 — U)g = imply 



Q(0)Q(P) 

1 (R- 1 
= 1 + - 



exp 

2 



R-l 
R+l 



B{U' - U) 



2\R+1 



p(\U'-U\ 2 ) § + 0(\R-l\ 3 ). 



Therefore, 



= l - 



1 f R-l 



2\R+1 



On the other hand, 
R-l 



exp 



R+l 



(3\U'-U\ 



P 2 (\u'-u\ 2 ) B 

+0(\R-1\ 3 ). (10) 



+\ (f^) 2 /? 2 (1^ - U\ 2 )- p + 0(\R- 1| 3 ). (11) 



Eq. JTUJ) and (JTTJ imply 



Q(/3) 2 



Q(P)QW) 



exp 



i?- 1 
i?+ 1 



l_|_±MGg) + 0(|i2-l| 3 ) 



and the proof is concluded. □ 
Proposition 1, while a powerful asymptotic result, has 
the disadvantage that it may produce negative num- 
bers for the acceptance probability in actual simulations. 
However, we can repair this in very straightforward fash- 
ion. Notice that the fact that Q(U) does not depend upon 
the temperature is not crucial for the proof of Proposi- 
tion 1. Thus, given any other well-behaved density of 
states Sl'[U, d((3)] depending perhaps on the inverse tem- 
perature through an adjustable parameter d(@) and such 
that 



&{\U-U% = M{p) 



(12) 



we still have 



R-l 



Ad '{13, f3') = 1 - -^—M (]3) + 0(\R~ 1| 3 ) . (13) 

In Eqs. (|12|) and (|13fl . the prime sign denotes the respec- 
tive averages or quantities obtained from their definition 
if Q(U) is replaced by Q.'[U,d{/3)]. From Eqs. © and 
(JT3J, we deduce 

Ac'(f3, f3') = Ac(f3, (3') + 0(\R- 1| 3 ) (14) 

for all il'[U,d{fi)\ that satisfy Eq. This sim- 

ple observation allows us to construct approximations 
Ac'(f3, f3') of order O (\R — 1| 3 ) for the acceptance prob- 
ability Ac(f3,f3') by considering a special functional de- 
pendence for the density of states W[U,d((3)], such that 
the resulting approximation is exact for a certain class of 
physical systems. 

We take this class to be the harmonic oscillators, for 
which Q'(U,d) = {2TT) d / 2 T{d/2)- l U d / 2 - 1 . We prove in 
Appendix I that for any d-dimensional harmonic oscilla- 
tor, 



Ac H {p,/3') = 



B(d/2,d/2) J a 



9 d/2-l( 1 _ e) d/2-l de 

(15) 



and 



M H (l3) 



-,2-d 



B(d/2,d/2)- 1 . 



(16) 



Here, B{d/2 1 d/2) denotes the respective Euler's beta 
function. 

For an harmonic oscillator, d is the dimension. For 
general systems, d = d{(3) is just a fitting parameter cho- 
sen such that Eq. $12\ is satisfied. In fact, Eq. (|12(1 . which 
in our case reads 

2 2 - d B(d/2,d/2)- 1 = M(f3), 

has always a unique solution because 2 2 ~ d B(d/2, cZ/2) — 1 
increases strictly from to +oo, as d also increases from 
to +oo. The following theorem is an immediate conse- 
quence of Eq. H14JI and of the discussion above. 

Theorem 1 (Incomplete beta function law) 

Consider an arbitrary thermodynamic system for which 
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M((3) is given as a function of temperature. Let d{(3) be 
the unique solution of the equation 



2 2 - d B(d/2,d/2)- L = M((3). 



Then, 
Ac(f3,f3') = 



1/(1+R) 



(17) 



B[d([3)/2,d(f3)/2] Jo 
x(l-6») !i(/5)/2 - 1 ^ + 0(|i?-l| 3 ), (18) 
with the error term canceling for harmonic oscillators. 



III. CONSEQUENCES OF THE INCOMPLETE 
BETA FUNCTION LAW 

There are several important consequences of Theo- 
rem 1. The first one concerns the development of an em- 
pirical law connecting the acceptance probabilities with 
the heat capacity and the ratio of the temperatures in- 
volved in the parallel tempering swap. The second one 
regards the optimal distribution of temperatures in par- 
allel tempering Monte Carlo simulations. Yet a third one 
is the statement that the efficiency of the swaps between 
neighboring temperatures decreases naturally as the in- 
verse square root of the dimensionality of the system. We 
analyze these consequences in some detail in the remain- 
der of the section. 



A. The empirical incomplete beta function law 

The property M(f3) is not usually determined in simu- 
lations, nor is it measured in experiments. It is therefore 
necessary to relate it to other thermodynamic properties, 
more precisely to the heat capacity. In addition, it is de- 
sirable to develop a version of the incomplete beta func- 
tion law involving the heat capacity rather then M(/3), 
even if the law has an empirical validity only. 

The quantity 

M(p)=0(\U'-U\) iJ 

measures the statistical fluctuations of the potential V(x) 
and is connected to the heat capacity of the system. More 
precisely, the Cauchy-Schwartz inequality says that 

32 /ITT' TT\\l ^ fl2 



However, 



p 2 (\U'-U\)j<(3 2 (\U' -U\ 



fF(\U'-U\')=J3 



1 



P " QW Jo Jo 



OO POO 



e -We-W 



n(U)Sl{U')(U' - UfdUdU' = 2f3 2 



QW Jo 



-pu 



xn(u)u 2 du 



- 0u Sl(U)UdU 



The last term in the equation above is twice the poten- 
tial contribution Cy (f3) to the total heat capacity of the 
system. (In this article, the heat capacity is always ex- 
pressed in units of the Boltzmann constant fee.) The 
total heat capacity sums both the potential and the ki- 
netic average square fluctuations and is given by the well- 
known formula 



C0) = C v (f3) + d/2. 
Then, the identity 

P 2 (\U'-Uf)_ = 2C v (f3) 

implies the inequality 

M(f3) 2 < 2C V ([3), 



(19) 



(20) 



(21) 



Eq. I|21|) suggests that M(j3) 2 is an extensive property 
of the physical system, property that may be very large 
for systems for which the heat capacity is also large. In 
fact, Sterling's formula shows that 

for large dimensional harmonic oscillators. The relation 
has a linear scaling with the dimensionality of the sys- 
tem (that is, with the number of particles). This scaling 
appears also for the heat capacity of harmonic oscillators 

C H , v 0)=d/2. 

For systems in condensed phase, for which an har- 
monic superposition is roughly a good approximation of 
the Boltzmann distribution, one may safely assume that 
the functional relationship between M((3) 2 and Cy(/3) is 
not very far from the one for the harmonic oscillator. Of 
course, this is always true in the low temperature limit. 
In this conditions, the solution d([3) of the equation 

2 2 - d B(d/2,d/2)- 1 = M{0). 

is approximately d(fj) = 2CV(/3). Replacing the last re- 
sult in Eq. I|18|) , we obtain the following empirical incom- 
plete beta function law: 



Ac((3,(3') 



1/(1+R) 



qC v (0)- 



B(C v (f3),C v ([3)) Jo 

x (l-0fv(/3)-i de _ (22 ) 



The good applicability of Eq. (122ft for realistic physical 
systems will be illustrated in Section IV for the case of 
a Lennard- Jones cluster. Eq. 1221) can be expressed in 
terms of the full heat capacity of the system with the 
help of Eq. I|19|l . The empirical incomplete beta function 
law is still exact for harmonic oscillators. 
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B. On the optimal schedule of temperatures for 
the parallel tempering simulation 

It is an empirical observation^ that the optimal sched- 
ule (i.e., the schedule for which all the acceptance proba- 
bilities for swaps between neighboring temperatures are 
equal) is given by a geometric progression of tempera- 
tures, if the heat capacity of the system is approxima- 
tively constant for the range \fi m in > Pmax] ■ The incom- 
plete beta function law gives a direct explanation of this 
phenomenon. As discussed in the previous section, the 
quantities M(f3) and Cv(f3) measure essentially the same 
information: average potential fluctuations and average 
potential square fluctuations, respectively. Thus, we ex- 
pect that M(/3) is also constant in regions where Cy(/3) 
is. As shown by Eq. I|18fl . m these conditions, the ac- 
ceptance probabilities are a function of the ratio of the 
temperatures involved in swaps only. Therefore, if the 
predetermined number of replicas is N (with N large 
enough so that Theorem 1 applies), then the optimal 
schedule is 



I < i < N, 



where 



R= (8 IB ■ W(w-i) 

Jt ' \t J max j Hmm J 



(23) 



(24) 



If one trusts it, the empirical incomplete beta func- 
tion law can be used to evaluate the optimal schedule 
for arbitrary systems in the following way. First, one 
computes a set of values for the potential part Cy (/3) of 
the heat capacity over the interval \fimin , Anax] using a 
geometric progression law. The results are then interpo- 
lated by cubic spline, for example. Only a rough estimate 
is necessary. Given a prescribed value p for the accep- 
tance probability and given the inverse temperature Pi, 
one computes /3i+i by solving the equation 



(l+ft+i/ft) 



B[C v {fl),Ov{fl)) Jo 

x (l-$) c vW- 1 d0=p, (25) 

where (3 — (fy + B i+ i)/2. One starts with f3\ — (3 m in and 
continues the procedure until the current inverse tem- 
perature becomes greater or equal to B max . This way, 
one determines both the optimal distribution of temper- 
atures and the minimal number of temperatures com- 
patible with the prescribed acceptance probability p. To 
ensure the validity of the approximation furnished by 
Theorem 1, the value of p should be large enough. In 
fact, values larger or equal to 0.5 are necessary anyway 
in order to have good mixing between the Monte Carlo 
walkers running at neighboring temperatures. 

Sure enough, one may use the full incomplete beta 
function law to find the optimal schedule. However, the 
computation of M(f3) requires extensive changes to the 
existing codes. In addition, as illustrated by the example 



described in Section IV, the empirical version of the in- 
complete beta function law may be sufficiently accurate 
for most systems of practical interest. The applicabil- 
ity of the law can also be tested during the computation 
of the heat capacity CV(/3), by comparing the observed 
values for the acceptance probabilities with the ones pre- 
dicted by the empirical incomplete beta function law. 



C. Loss of efficiency with increasing dimension for 
the parallel tempering method 

In this subsection, we show that the minimum num- 
ber of intermediate temperatures N(d,p) that ensures 
an acceptance probability greater or equal to some preset 
value p € (0, 1) for a d-dimensional system increases nat- 
urally at least as the square root of the dimensionality for 
condensed-phase systems (both solids and liquids). This 
observation was first made by Hukushima and Ncmoto. 20 

We begin with a rigorous mathematical analysis for 
harmonic oscillators. For them, the optimal schedule is 
a geometric progression [because M{B) is independent of 
temperature] and the minimum value of N(d,p) is given 
by 



N H (d,p) 



\n(B 

max/ ftmin) 



HR(d, P )} 



(26) 



where [x] is the integer part of x and R(d,p) is the solu- 
tion of the equation 



0d/2-l^_ e y/2-l d6 = p 



B(d/2,d/2) Jo 

As shown by Theorem 2 of Appendix 2, 

i 

2 f ! + «(■*. p) 

d^o B{d/2,d/2) J 



12-1 



(l-0)d/2-l d6 =p 



if and only if 
lim Vd 



1 + R(d,p) 



where erf -1 is the inverse of the error function. Straight- 
forward manipulations show that the last relation is 
equivalent to 

R(d, p) = 1 + ^lerF 1 (1 -p)+o (d- 1 ' 2 ) , 



or to 



HR{d,p)] 



Vd 
2V2 



erf 



\l-p)+o(d- 1 / 2 ). (27) 



We remind the reader that, in general, the notation Xd — 
x + o{d a ) means 

lim (x d - x)/d a = 0. 
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From Eqs. (JSJjJ and l(2*7|) . we learn that 

N H (d, P ) =d l/2 ^H^aJf3 mm ) +(} ( d l/2\ (28) 

4crF (1-p) V / 

relation that is similar to Eq. (9) of Ref. ITsl 

For condensed-phase systems, the heat capacity is usu- 
ally larger than what an harmonic approximation pre- 
dicts and, therefore, so is M(f3). Consequently, for such 
systems, the acceptance ratio for harmonic oscillators 
Ach(P,/3') is usually not attained. In this conditions, 
all we can say is that N(d,p) must satisfy the relation 



N.„ 



13 is, of course, the number of particles in the 



N(d,p) > d 1/2 



v2 \n(j3max I ftmi'i 

4erf~ 1 (l -p) 



(29) 



inequality which demonstrates our claim that there is a 
curse of dimension for the parallel tempering method. 
In this respect, notice that running N{d,p) independent 
replica requires N(d,p) times more computational power. 
However, given the improvement in the quality of the 
sampling brought in by parallel tempering, this loss of 
efficiency is an acceptable price to pay. 



IV. A NUMERICAL EXAMPLE 

In this section, we verify the validity of the empiri- 
cal incomplete beta function law for the neon realiza- 
tion of the LJ13 cluster. This example is representative 
of the type of applications one is likely to encounter in 
practice. We shall also illustrate the use of the incom- 
plete beta function law for the design of optimal sched- 
ules for parallel tempering. With the help of the iden- 
tity C(/3) = CV(/3) + d/2, the empirical incomplete beta 
function law is expressed as a functional of the full heat 
capacity as 

Ac(/3, /?') » 2 B(C{P) - d/2, C(/3) - d/2) _1 
r l/(l+R) 

X / d C(0)-d/2-l( 1 _^C(0)-d/2-l d6i (30) 

JO 

where d — 3 • 13 = 39 is the dimensionality of the cluster. 
The total potential energy of the Nei3 cluster is given 

by 



N p N p 
i<j 



(31) 



i=l 



where Vtjfoj) is the pair interaction of Lennard- Jones 
potential 



V L j{r lj ) = 4:6 l, j 



and V c (i"i) is the confining potential 



12 / x 6 
0~LJ \ / LJ 



V c {n) = e LJ 



|r t - R c 



211 



(32) 



(33) 



system. The values of the Lennard- Jones parameters o~lj 
and clj used are 2.749 A and 35.6 K, respectively. The 
mass of the Ne atom was set to mo = 20, the rounded 
atomic mass of the most abundant isotope. R cm is the 
coordinate of the center of mass of the cluster and is given 
by 



Rc 



N p 



(34) 



Finally, R c — 4a lj is the confining radius. The role of 
the confining potential V c (i"i) is to prevent atoms from 
permanently leaving the cluster since the cluster in vac- 
uum at any finite temperature is metastable with respect 
to evaporation. 

The range of temperatures employed was 3 K to 30 K. 
In this range of temperatures, the system undergoes two 
phase transitions: a solid-liquid transition at about 10 K 
and a liquid-gas transition at about 18 K, respectively. 
Surely, the terminology we employ is more or less abuse of 
language, because true first-order phase transitions hap- 
pen only in the limit of an infinite number of particles. 
However, Fig. ^ clearly shows two pronounced maxima 
in the heat capacity of the system, maxima that separate 
the three phases. In the limit of infinite number of parti- 
cles, these maxima sharpen to a delta function and their 
temperature value is usually lowered to the correspond- 
ing bulk values. 

The parallel tempering Monte Carlo simulations are 
carried out using a total of 32 parallel streams, each 
running a replica of the system at a different temper- 
ature. The temperatures range from 3 K to 30 K and are 
distributed in geometric progression. For each stream, 
the coordinate sampling is performed with the help of 
the Metropolis algorithm^, The basic Monte Carlo steps 
consist in attempted moves of the physical coordinates 
associated with a given particle. Each attempted move 
is then accepted or rejected according to the Metropolis 
logic. By attempting to move the particles one at a time, 
we can increase the maximum displacements and ensure 
better quality of sampling. The maximum displacements 
are adjusted in the equilibration phase of the computa- 
tion, so that each of the acceptance ratios eventually lies 
between 40% and 60%. The 32 (statistically) indepen- 
dent streams of random numbers necessary for the simu- 
lation are obtained with the help of the Dynamic Creator 
package 

As a counting device, we define a pass as the mini- 
mal set of Monte Carlo attempts over all particles in the 
system. Because we update the neon atoms in successive 
fashion, a pass consists of 13 basic steps. One also defines 
a block as a computational unit made up of 100 thousand 
passes. The size of the block is sufficiently large that the 
block averages of the various quantities computed are in- 
dependent for all practical purposes. The accumulation 
phase of the simulation has consisted of 100 blocks for a 
total of 10 million passes per temperature. The accumu- 
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FIG. 1: Heat capacities (in units of fee) as a function of the 
temperature T (in Kelvin). The data for the plot have been 
computed with the help of the optimal schedule determined 
in Section IV. The error bars (twice the standard deviation) 
are smaller than the plotting symbols. 



FIG. 2: Observed acceptance ratios as a function of the tem- 
perature T (in Kelvin). A geometric progression schedule has 
been utilized. The plot is in direct relationship with the heat 
capacity curve (see Fig. 0, as predicted by the incomplete 
beta function law. The error bars are about the size of the 
plotting symbols. 



lation phase has been preceded by an equilibration phase 
of 20 blocks. 

Parallel tempering swaps between neighboring temper- 
atures are attempted every 100 passes in an alternat- 
ing fashion (first, with the closest lower temperature and 
then, with the closest higher temperature). The only 
exceptions are the two end temperatures, which are in- 
volved in swaps every 200 passes, only. The acceptance 
probability of swaps aci at a given temperature fli is com- 
puted as the fraction of accepted swaps involving that 
temperature. Thus, except for the end temperatures, the 
computed values are equal to 

ac, = - [Ac(ft_i,A)+^c(ft,A+i)] ■ 

Because the intermediate temperatures are distributed 
in geometric progression, R = (30/3) 1 / 31 is a constant. 
In this case, the empirical incomplete beta function law 
given by Eq. (|5U|l says that the dependence of the accep- 
tance probabilities with the temperature is a functional 
of the heat capacity only. This observation is well sup- 
ported by the observed acceptance probabilities, which 
arc plotted in Fig. 

During the Monte Carlo simulation, besides acceptance 
ratios, we have also evaluated the heat capacities at dif- 
ferent temperatures. Then, the heat capacities have been 
interpolated using a cubic spline and acceptance ratios 
ac^ have been computed with the help of the empirical 
incomplete beta function law. More precisely, we have 



computed 

ac^ = i [Ac' + Ac'{(3 h (3 i+1 )] , 

where Ac'(j3 : j3') is given by the right-hand side of 
Eq. (|30fl . In Fig. [31 we plot the absolute values of the 
differences between the observed and computed accep- 
tance probabilities. These values are smaller than the 
estimated error bars for the same differences. In fact, 
the maximum difference between the two curves is less 
than 0.008 and therefore, for all practical purposes, the 
empirical version of the incomplete beta function law is 
correct. 

As discussed in Section III.C, the heat capacity for 
solids is usually larger than what the harmonic ap- 
proximation predicts, except for the low temperature 
limit. The liquid phase has an even larger heat capacity, 
whereas the gas phase has a smaller one. In the high- 
temperature limit, the Cy(/3) component of the heat ca- 
pacity decreases to zero. As a consequence, the accep- 
tance probabilities increase to 1. At the critical points, 
the acceptance probabilities for parallel tempering swaps 
decrease very much, as the heat capacities increase sud- 
denly. This analysis is consistent with the results shown 
in Fig. |21 Therefore, the loss of efficiency for solids and 
liquids is even more pronounced than the dT x l' 2 decay 
computed for d-dimensional harmonic oscillators. 

Using the strategy described in Section III.B, we have 
determined the optimal schedule of temperatures nec- 
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FIG. 3: The dot line shows the absolute value of the differ- 
ences between the observed values for the acceptance prob- 
ability of swaps and the ones predicted by the empirical in- 
complete beta function law. Also plotted (dash line) are the 
estimated error bars for the differences \aci — ac'A. 



essary to achieve a constant acceptance probability of 
p = 0.75 for all swaps between neighboring temperatures. 
The minimal number of temperatures needed has been 
found to be 34 (only the temperatures located in the in- 
terval [3 K, 30 K] are counted). Then, we have performed 
a second Monte Carlo simulation to verify the validity of 
the schedule. The plot in Fig. 01 demonstrates that the 
computed schedule works very well. This explicit appli- 
cation illustrates the utility of the empirical incomplete 
beta function law for the determination of optimal tem- 
perature schedules in parallel tempering simulations. 



V. SUMMARY AND CONCLUSIONS 

We have successfully and rigorously related the ac- 
ceptance probabilities for parallel tempering swaps to 
the ratio between the temperatures involved in the swap 
and the average statistical fluctuations of the potential 
at some intermediate temperature. The respective law, 
called the incomplete beta function law, is exact for 
harmonic oscillators and of order 0(\(3' — /3| 3 ) for arbi- 
trary systems. We have also demonstrated that there is 
a loss of efficiency in parallel tempering simulations of 
condensed-phase systems with the increase of dimension- 
ality. The loss of efficiency is at least d -1 / 2 , the value 
computed for harmonic oscillators. 

Motivated by the fact that the existent Monte Carlo 
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FIG. 4: Observed acceptance ratios for the optimal schedule 
of temperatures determined with the help of the empirical 
incomplete beta function law. The deviations from the ideal 
result of p = 0.75 are comparable to the statistical noise. 



codes do not allow for the computation of the average 
potential fluctuations without extensive reprogramming, 
we have developed and tested the empirical incomplete 
beta function law. This empirical law connects the accep- 
tance probabilities of the parallel tempering swaps with 
the heat capacity of the system. The law has been exten- 
sively verified for the LJ13 cluster, on a range of temper- 
atures that spanned three thermodynamic phases. The 
empirical incomplete beta function law provides a direct 
justification of the observation that a geometric progres- 
sion is the optimal schedule for systems and regions of 
temperatures where the heat capacity is almost constant. 
Finally, the use of the empirical incomplete beta function 
law for the construction of optimal temperature schedules 
has been demonstrated. 

As opposed to its empirical version, the incomplete 
beta function law given by Theorem 1 is an exact math- 
ematical statement, valid for all systems asymptotically, 
for close enough temperatures. For strongly anharmonic 
systems (for instance, systems for which the sampling 
is performed on a lattice, such as spin glasses and self- 
avoiding random walks), the empirical version of the in- 
complete beta function law may fail. In such cases, the 
rigorous incomplete beta function law should be used for 
the development of optimal temperature schedules. As 
discussed in the text, the evaluation of the average 
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l d 

x/3|V(x')-y(x)| (35) 



9 



requires however a double integral over the configuration 
space. This integral can be computed by doubling the 
number of temperatures in the parallel tempering sched- 
ule, as follows 



00 = 01 < 02 = 03 < 04 = 05 < ■ ■ ■ < 02N-2 = 



2N-1- 



Then, to compute M(/?j), one collects the values of 
the differences |V(x) — V(x')| any time a swap be- 
tween equal neighboring temperatures 0i and 0%±i is at- 
tempted. Of course, swaps between equal temperatures 
are always accepted. The values M(0i) are interpolated 
by a cubic spline. The determination of the optimal tem- 
perature schedule then proceeds in a way similar to the 
approach utilized in Section IV. 

While the reader may object that the introduction of 
an intercalated set of identical temperatures is an ad- 
ditional computational burden, many times, there are 
certain advantages in doing so. For large dimensional 
systems, coupled independent replica running at identi- 
cal temperatures constitute an elegant device for paral- 
lelizing the Monte Carlo code, whenever the number of 
available compute nodes is at least twice the number of 
temperatures in the parallel tempering schedule. Nowa- 
days, with the advent of inexpensive cluster computing, 
this is the case with many research groups. Furthermore, 
in this setting, the computation of the heat capacity and 
other properties of the system can be done by means of 
unbiased estimators, as shown by Eq. l]2Up. 

On a more general level, we hope that a better under- 
standing of the laws governing the acceptance probabili- 
ties for swaps in parallel tempering methods may lead to 
useful research in improving the efficiency of the meth- 
ods. In the meantime, we recommend that the dimen- 
sionality of the systems be maintained as low as possible, 
for instance, by adiabatically reducing those degrees of 
freedom that do not lead to significant degradation in the 
quality of the final results. 
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APPENDIX A: EVALUATION OF INTEGRALS 

For an arbitrary d-dimensional harmonic oscillator 
whose global minimum is zero, the density of states is 
given by the formula Q(U) = (2-K) d ' 2 U d / 2 - x /Y(d/2). 
In these conditions, it is but a simple exercise to show 



that the acceptance probability for the parallel temper- 
ing swaps is given by the equation 

*c H[ p,0) T(d/2)2 J Q J q 

xU'^min^eV'-Mv'-v^dUdU'. (Al) 

We also want to evaluate the quantity Mh(0) [see 
Eq. 116|) ]. which is given by the formula 

ad+1 rco poo 

xU ,d/2-ip> -U\dUdU'. (A2) 

While for an harmonic oscillator the parameter d is an 
integer, we compute the two integrals above under the 
assumption that d is a strictly positive real number. We 
prove that the value of Ach(0, 0') is given by the incom- 
plete beta function 



Ac H (0,0') = 



1/(1+*) 



B(d/2,d/2) J 
where B(a, b) is Euler's beta function 



9 d/2-l (1 _ g yl/2-l d ^ 

(A3) 



B(a,b) = f a - 1 (l-, 
Jo 



\b-l 



d6 



and 



R = max 



i. 1 

/?' 0' 



In addition, we prove that 

M H {0) = 2 2 - d B(d/2,d/2)-\ 



(A4) 



Because the function Ach(0,0') is symmetrical in its 
arguments, we may assume without loss of generality that 
0' > 0, so that R — 0'/0>\. By decomposing the do- 
main of the integral against U' in Eq. (|A1|I in two regions 
with U 1 < U and U' > U respectively, it is straightfor- 
ward to see that 



Ac H (0, 0') = l- 1(0, 0') + 1(0', 0), (A5) 



where 



1(0,0') = 



Ad/2 poo r U 



dU / dU'e^ u 



(00') d / 2 

T(d/2) 2 J - J 

xe -f3'U' u d/2-l u/d /2- K (A6) 



Performing the substitution U' = UQ, we obtain 

od r?d/2 poo pi 
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Performing the change of variables 9 = R x (t 1 — 1), we 
conclude 



1 



B(d/2,d/2) J 1/(R+1 



B{d/2,d/2) Jo 



t d/2-i( 1 _ t y/2-i dt 

t d/2-l( 1 _qd/2-l d ^ ( A8 ) 



where 



B(d/2, d/2) = = f t d/2 ~\l - t) dj2 - x dt 

is Euler's beta function. The value of I(/3',/3) is ob- 
tained by replacing R with 1/R in the first expression 
of Eq. (|A8|) . We compute 



/(/?',/?) 



1/(1+R) 



t d/2-l {l _ t) d/2-l dt 



B(d/2,d/2) y 

Eq. 1A.3|) follows easily from Eqs. l|A5jl . l|A8jl . and l|A9jl . 

after easy simplifications. 

If applied to Eq. (| A2|) . the same decomposition and 
coordinate transformations used in the proof of Eq. (|A3|) 
lead to the equation 

M H {0) = 2 r(d ^^ / (26 - l)0 d / 2 - x (l - 9) d / 2 - l d9 



T{d/2f 



1/2 



2dB(d/2,d/2)~ 1 



2I 1 --B(d/2,d/2) 



(A9) 



where 



h = / o d ' 2 {i - ef/^de. 

Jl/2 

Integrating by parts the last integral, we obtain 



h 



d2 d 



+ 0d/2-i(i-e) d / 2 de = — + i 2 , (Aio) 



1/2 



d2 d 



where 



1/2 



(l-0)d/2-l 6 d/2 de _ 



Combining Eq. IjAlOfl with the identity 

h+h= B(d/2, d/2 + 1) = ^B(d/2, d/2), 



we obtain 



2h 



-B(d/2,d/2), 



d2 d ' 2 

which, after replacement in Eq. I)A9|1 . produces Eq. 1A4|I . 
APPENDIX B: A LIMITING THEOREM 

Theorem 2 Let {ctd}d>i be a sequence of positive num- 
bers convergent to a > 0. Then, 



/2-l (1 _0)d/2-l d 

= 1 - erf(2 1 / 2 a), (Bl) 



lim , 

a— »oo B(d/2,d/2) J 



where 

erf(x) 

is the error function. 
Proof. We compute 
2 



e~* dt 



= 1 - 



= 1 - 



B(d/2,d/2) J Q 
2 



B(d/2,d/2)ji_j^ 



2 d - 2 B(d/2,d/2) J 



9 d/2-l( 1 _ 0) d/2-l d6 
d/2-l( 1 _ 0) d/2-l d9 

d (1 - At 2 ) d / 2 - l dt, 



where we have used the transformation of coordinates 
9 = 1/2 — t for the last expression. Next, we perform the 
substitution t = Odd/d 1 / 2 and obtain 



I r , = l 



4a„ 



d 1 / 2 2 d ~ 1 B(d/2,d/2) J 
Sterling's formula implies 

1 



2„,2\ d/2-1 



A6 2 at 



d6 



J^o B(d/2,d/2)2 d - 1 d 1 / 2 

r(d) i i 



lim 



d^ooT{d/2) 2 2 d - 1 d 1 / 2 V2T 
Therefore, remembering that a n — > a, we obtain 

Aa f 1 ( 29 2 a 2 \ d/2 ~ 1 

lim I d = 1 - -= lim / 1 - — ^ 1 d# 

rf^oo v27T d— too Jo \ U/Z 



The equality 



d/2 



the fact that the above sequences are bounded by 1 for 
all 9 <E [0,1], and the dominated convergence theorem 
imply 



lim 

d — too 



Ofl2„2\ d/2-1 r i 
20 a d\ M _ f ^-2e 2 a 2 d ^ 



d/2 



d9= I r 
'o 



Thus, 



lim I d = 1 = ; < 

d^oc ^/2n Jo 



4« / - 2a ^ df) 



1- / Q e -* 2 ( it = l-erf(2 1 /2 a ) 



and the theorem is proven. 
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